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Abstract 

We consider a simplified coarse-grained model for colloid-polymer mixtures, in which polymers 
are represented as monoatomic molecules interacting by means of pair potentials. We use it to 
study polymer-colloid segregation in the presence of a quenched matrix of colloidal hard spheres. 
We fix the polymer-to-colloid size ratio to 0.8 and consider matrices such that the fraction / of the 
volume that is not accessible to the colloids due to the matrix is equal to 40%. As in the Asakura- 
Oosawa-Vrij (AOV) case, we find that binodal curves in the polymer and colloid volume-fraction 
plane have a small dependence on disorder. As for the position of the critical point, the behavior is 
different from that observed in the AOV case: while the critical colloid volume fraction is essentially 
the same in the bulk and in the presence of the matrix, the polymer volume fraction at criticality 
increases as / increases. At variance with the AOV case, no capillary colloid condensation or 
evaporation is generically observed. 

PACS numbers: 61.25.Hq, 82.35.Lr 



I. INTRODUCTION 



The study of the properties of fluids in porous materials jl| is very important for its 
technological applications in many fields of science, from physics and chemistry to geology, 
engineering, and even agriculture. In this paper we consider the demixing of a binary mixture 
of nonadsorbing colloids of radius Rc and of polymers of radius of gyration Rg, with the 
purpose of understanding how the porous structure influences the phase behavior. We are 
considering mesoporous disordered materials, like silica gels, which present large pores and 
can thus adsorb mesoscopic particles like colloids. In the bulk the phase behavior of colloid- 
polymer mixtures depends on the parameter q = Rg/R^.. For q small, q < qc {qc ~ 0.2-0.3 
under good-solvent conditions), only a fluid-solid transition is observed, analogous to that 
observed in hard-sphere systems. For q > qc also a fluid-fluid transition occurs, between a 
colloid- liquid (polymer-poor) and a colloid- gas (polymer-rich) phase. It is such a fluid-fluid 
transition that will be investigated in the present paper. 

At least qualitatively, many aspects of the behavior of polymer-colloid mixtures can be 
understood by using the Asakura-Oosawa-Vrij (AOV) model , which gives a coarse- 
grained (CG) description of the mixture. Polymers are treated as an ideal gas of point 
particles, whose radius is identified with the radius of gyration Rg, which interact with 
the colloids (hard spheres of radius Rc) by means of a simple hard-core potential. This 
model is extremely crude since it ignores the polymeric structure and polymer-polymer 
repulsion, which is relevant in the good-solvent regime. Nonetheless, it correctly predicts 
polymer-colloid demixing as a result of the entropy-driven effective attraction (depletion 



interaction) between colloidal pairs due to the presence of the polymers |4ll2|. AOV colloid- 
polymer mixtures in a porous matrix have been studied in Refs. |l3l-ll8| by means of density- 
functional theory, integral equations, and Monte Carlo (MC) simulations. The nature of 



the critical transition has been fully clarified 



14j-|l7l|: if the obstacles are random and there 



is a preferred affinity of the quenched obstacles to one of the phases, the transition is in 
the same universality class as that occurring in the random-field Ising model, in agreement 
with a general argument by de Gennes 19|. If these conditions are not satisfied, standard 



Ising or randomly dilute Ising behavior is observed instead, see Refs. 
considered the AOV model and investigated 
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2l|. Recently, we 



18| how demixing is infiuenced by the amount 



of disorder and by its nature. We found that demixing was, to a large extent, only dependent 
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on the fraction / of the volume that is not accessible to the colloids due to the presence of the 
random matrix. The matrix topology was instead largely irrelevant. Moreover, we observed 
the possibility of capillary condensation of the colloids: for some values of the parameters, a 
colloid-gas bulk phase is in equilibrium with a colloid-liquid phase in the matrix. 

The AOV model completely neglects polymer-polymer interactions, hence it may only 
be quantitatively predictive close to the 9 point where, to some extent, polymers behave as 
ideal particles. In this paper we make a first step towards the inclusion of polymer-polymer 
interactions, allowing us to study systems under good-solvent conditions. We still use a CG 
model in which polymers are treated as monoatomic molecules, but we include a polymer- 
polymer repulsive pair potential which is such to reproduce the correct thermodynamics in 
;he low-density limit. In recent years, this class of CG models has been extensively studied 
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24j . It is now clear that, unless one includes many-body interactions or considers density- 



dependent potentials, they are quantitatively predictive only in the dilute regime in which 
the polymer volume fraction rjp [rjp = Cp/c*, where Cp = Np/V is the concentration and 
c* = 3/(47ri?^)] satisfies rjp < 1. If g < 1 (the so-called colloid regime), fiuid-fiuid demixing 
occurs for f]p ^ 1, hence monoatomic CG models are expected to provide reasonably accurate 
results. If, instead, g > 1 (the protein regime), demixing occurs for rjp > 1, hence one must 
use more sophisticated CG multiblob models 25|] in which each polymer is represented by 
a polyatomic molecule and each atom corresponds to a polymer blob of size of the order 
of that of the colloid. In this paper, we investigate the case q = 0.8, as in our previous 



work 



1^ 



. Since colloids are larger than polymers, a CG approach based on a single-blob 



representation of the polymers should be adequate. 

The exact polymer-polymer and polymer-colloid pair potentia 
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26 



s appropriate for single- 



-|29|. The polymer-polymer 



blob models have been computed in several papers {4, 
potential has a Gaussian shape with a width of the order of Rg and is significantly different 
from zero {upp{r) / Upp{0) > 0.01) up to r ^ 2.5Rg. Analogously, the typical colloid-polymer 
potential has a tail, which is still significant for r ^ '^{Rg + Rc)- The presence of these 
tails makes simulations quite slow, since one must consider a large number of neighboring 
molecules in each updating step. Since our simulations are extremely complex and time- 
consuming — we must average over a large number of disorder realizations — we decided to 
replace the exact potentials with simpler ones that generalize the AOV interactions. They 
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are square potentials, hence have no tails, and thus allow us to determinate the energy quite 
fast. Of course, this simplification implies that our results cannot be quantitatively accurate. 
Still, our simulations should allow us to understand how the presence of the porous matrix 
changes the behavior of the polymer-colloid mixture under good-solvent conditions. 

The paper is organized as follows. In Sec. [ITlwe present the simple model which we will use 
and give the basic definitions. In Sec. IIIII we discuss the model in the bulk, while in Sec. HV] 
we determine the binodal lines and the critical-point positions for polymers and colloids 
adsorbed in two different colloidal matrices. Finally, in Sec. |V] we present our conclusions. 



II. DEFINITIONS 



A. Models 



In this paper we model polymers as soft effective spheres. Polymers of radius of gyration 
Rg are represented by monoatomic molecules interacting with pair potential 

{epp for r < aRg, 

(1) 
for r > aRg, 

where a and epp are parameters which are fixed below. Since we expect demixing to occur 
for rjp < 1, we have fixed a and e^p to reproduce accurately the compressibility factor 
Z = p/{kBTc) {p is the pressure and c the concentration) in this density interval. In 
practice, we have determined the second virial coefficient -B2,pp and Z at r]p = 1 using model 
([T]) for several values of the parameters and we have compared the results with those obtained 



in full-monomer simulations 
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31| . Requiring the model to reproduce the estimate [30 1 



^2,pp = B2^ppRg ~ 5.50 and the estimate |31f] of Z for r^p = 1, we obtain 

a = 1.58 epp = 1.096. (2) 

It is important to note that, although the thermodynamics is quite well reproduced up to 
rjp ^ 1 and with small errors up to ?7p = 2 (at such value of rip the difference between Z 
computed in the present model and that for polymers is 8%), the intermolecular structure 
is poorly reproduced. For instance, in this model the intermolecular distribution function 
is discontinuous and oscillates due to the discontinuity of the potential, a behavior which is 
not observed in polymer systems. 



Let us now introduce the colloids. Two colloids interact with a hard-sphere potential 

for r < 2Rr 




Uccir) = { (3) 

for r > 2Rc, 

while colloid-polymer interactions are modelled by taking a soft version of the AOV hard-core 
potential 

/ N ) for r <Rc + Eg, 

Ucp{r) = { (4) 

for r >Rc + Eg. 

The parameter ecp{q) is fixed so that the thermodynamics is exactly reproduced in the 
low-density limit. For this purpose we compute the universal polymer-colloid second- virial 
combination 

A,,., = B,,^,Ef/'Ef/' = lE-'/'Ef/' J dh (1 - e-«-('-)) = 

= f (l-e--^)(l + g)V'/^. (5) 

Estimates of A2,cp were obtained in Ref. [29i] for polymers under good-solvent conditions. 
We fix Scpiq) so that the value of A2^cp in our model is the same as that obtained in Ref. 



Results are reported in Table [H As expected, ecp(g) diverges as g — )■ 0, since in this limit 
polymers are quite well described by hard spheres. In the opposite limit instead, the inter- 
action energy vanishes, a phenomenon which is related to the fact that, for large values of 
q, the polymer can wrap around the hard spheres. Hence, in the CG model in which each 
polymer is replaced by a monoatomic molecule positioned in the polymer center of mass, 
there is a significant probability that the CG polymer and the colloid are one on top of 



the other. If we compare potentials (jlj) with the exact ones reported in Refs. |27l . |29| . we 
observe that Ucp{r), which is essentially an average potential, overestimates the interaction 
energy for Eg < r < Ec + Eg, while it significantly underestimates Ucp{r) close to overlap 
(for instance, the correct Ucp{r) diverges for r — )■ if g < 1). Note that a more accurate 
model could have been obtained by also introducing a parameter acp to specify the range of 
the polymer-colloid interactions, as in Eq. ([1]). However, to fix an additional parameter we 
would have needed some additional thermodynamic information (for instance, the pressure 
for a finite value of the polymer and colloid densities), which was not available to us. 
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5 0.346 
2.5 0.725 
1.25 1.363 

1 1.642 
0.8 2.035 

TABLE I: Estimates of the parameter e, 



cp- 



We mention that a different simphfied model including polymer-polymer interactions was 
introduced in Ref. [32| and studied for q = 0.8. However, with their parameter choices, 
thermodynamics is not reproduced in the low-density limit. Indeed, if we consider the 
second virial coefficients, their model gives A2,pp ~ 1.79 and A2^cp ~ 17.2 to be compared 
with estimates A2,pp ~ 5.50 and A2^cp ~ 14.8 obtained from direct polymer simulations 
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30|. While the colloid-polymer coefficient is reasonably close to the correct one (they 



differ by 15%), the second virial coefficient for polymers is significantly smaller. Hence, the 



model of Ref. 32| appears to be more appropriate to describe polymers in the crossover 
region close to the 6 point than in the good-solvent regime. 

The colloidal matrix has been introduced as in our previous work. We consider a random 
distribution of quenched hard spheres of radius -Rdis- The matrix-colloid and matrix-polymer 
interaction potentials are therefore: 



Ucd[r 




for r < i?dis + Rc, 
for r > Rdis + Rc, 



for r < Rdis + Rg, 
Upd{r) = \ (6) 

for r > i?dis + Rg, 

where gdis = Rg/Rdis- For _Rdis 0, i.e., for gdis — ^ oo, we have ecq{qdis) 0. As we already 
discussed, this reflects the fact that the polymer can easily wrap around the small quenched 
colloid, which implies that the CG polymer can easily overlap with it. As a consequence of 
this, the matrix becomes less and less repulsive as gdis increases, so that in the limit i?dis — ^ 0, 
i.e. gdis oo, we obtain bulk behavior. This is, of course, an artifact of the model, which 
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3, 



in which a single polymer is 



can only be eliminated by using multiblob approaches 
represented by many blobs, whose size is of the order of that of the quenched colloid. 

In the simple model we consider, disorder is characterized by two parameters, the reduced 
concentration c = Cdis-Rc (^dis = ^dis/^, where A^dis is the number of quenched hard spheres 
present in the volume V) and the ratio -Rdis/ Rc- However, as we already discussed in Ref. 



it is much more useful to characterize the amount of disorder by using the volume fraction 
/ which is not accessible to the colloids due to the matrix. To define it precisely, consider 
the region TZ in which the (centers of the) colloids are allowed: 

n= {r ■.\v-Ti\> R, + i?dis, for all 1 < i < A^dis }, (7) 

where rj is the position of the i-th hard sphere belonging to the matrix. If Vn is the volume 
of the region 71, we define 



V ' 

where [Vji] is the average of Vn over the different matrix realizations. 
B. Simulation details 

In this work we investigate the effect of disorder on the fluid-fluid binodals for q = 0.8 



which is the case we investigated in our previous work 18|]. We perform simulations in the 
absence of the porous matrix and for / = 0.4, i?dis/-Rc = 0.2, 1.0. In order to determine the 
coexistence curves we combine the grand-canonical algorithm with the umbrella sampling 



33| and the simulated-tempering method [SJ], as discussed in the Appendix of Ref. |18 |. 
The grand partition sum for each disorder realization is 

EiV, zp, z,) = J2 ^p'^c'QiV, Np, AT,), (9) 

where Q{V, Np, Nc) is the configurational partition function of a system of A'p polymers and 
Nc colloids in a volume V, and Zp and Zc are the corresponding fugacities. In Eq. ([9]) we 
normalize QiV, Np, Nc) so that QiV, 1, 0) = QiV, 0, 1) = V, hence Zp and Zc are dimensionful 
parameters. 

The system shows, both in the bulk and in the matrix, a demixing transition. For Zp < 
^p,crit a single phase exists, while for Zp > Zp^crit coexistence is observed along the hne Zc = 



8 



z*{zp). In the MC simulations the position of the demixing curve can be determined by 
studying the disorder averaged histograms of Nc and Np, which are defined as 



(10) 
is the 



where 6{x,y) is Kronecker's delta [6{x,x) = 1, S{x,y) = for x 7^ y], {■)gczpZc 
grand-canonical ensemble average, and [■] is the average over the matrix realizations. In 
the two-phase region the histograms show a double-peak structure. In order to obtain z* at 
fixed Zp in a finite volume, several different methods can be used. We followed two different 



recipes, the equal-area and the equal- height methods, as discussed in Ref. [l8|. They give 
completely consistent results: the results we report have been obtained by using the equal- 
height method. 



III. BULK BEHAVIOR 



Before considering the model in the presence of the matrix, we determine the phase 
behavior in the bulk. We use the algorithm described in Ref. ISj. One Monte Carlo iteration 
consists in 3 simulated-tempering fugacity swaps and 1000-5000 grand-canonical moves in 
which colloids and polymers are inserted or removed. For each value of Zp, we perform 
A^ini iterations to determine the umbrella functions and then Alitor iterations to measure the 
histograms. Typically, A^ini varies between 5000 and 20000 N^, while Niter is of the order 
of 5 ■ 10^ Nm- Here is the number of colloid fugacities which are sampled together in the 
simulated-tempering simulation; we take Nm ~ 10. 

The demixing curves have been determined for L/Rc = 14 and L/Rc = 16, to identify 
size effects. We have also performed simulations for L/Rc = 13 close to the critical point, to 
better determine its position. In Fig. [1] we report the liquid-gas coexistence curve in terms 
of the colloid and polymer volume fractions rjc and rjp defined as 

4 



4 „s 
Vc = ^-^CcR^, 



Vp = ■^'^CpRg, 



(11) 



where Cc and Cp are the concentrations of the colloids and of the polymers, respectively. It 
is clear that size effects are quite small, indicating that our results provide a good estimate 
of the true (infinite- volume) binodal curve. As already observed in Ref. [4|, the presence 
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FIG. 1: (Color online) Bulk demixing curve (solid symbols) and coexistence diameter (empty 
symbols) in the (r/c^/p) system representation for L/Rc = 14,16. We also report the demixing 



curve of the AOV model, solid line, from Ref. 
point. 



18l |: the solid circle gives the position of the critical 



of polymer-polymer interactions significantly changes the binodal curve: the values of rjp at 
which demixing occurs are significantly larger in the presence of polymer interactions than in 
the AOV case. An unusual feature of our results for the coexistence curve is that rjp along the 
coexistence line in the colloid-liquid phase slightly increases as rjc increases beyond 0.35. This 
feature is not observed in the more accurate model used in Ref. {4], nor in experiments. It is 
probably an artifact of the model, a consequence of the discontinuous nature of the potentials 
or of the fact that polymers and colloids can overlap paying a relatively small energy penalty 
of the order of 2kBT, a phenomenon which is not possible in the exact representation. 

We also report the liquid-gas demixing curve in the reservoir representation. We report 
the results in terms of the (ideal) reservoir packing fraction 

r^7^ ^ ^ZpRl, (12) 

which is the quantity used in the AOV model, and in terms of the reservoir packing fraction 
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FIG. 2: (Color online) Bulk demixing curve in the reservoir representation for L/Rc = 14, 16. r]p' 
is the ideal reservoir packing fraction, while r]p is the true reservoir packing fraction. 



r,id 



Tjp, which is the volume fraction of the polymers in the absence of colloids at a given Zp. We 
have determined it by grand-canonical simulations setting = 0. The two quantities r/^'^'^ 
and Tjp are related by 

v;'"' = v;exp[Pf,^;^^\r^;)], (13) 

where ^Jip'^^\rip) is the polymer excess chemical potential which can be computed by using 
the equation of state. In Fig. [2] we report the liquid-gas coexistence curve by using both 
definitions. Size effects are small also in this case, except close to the critical point. It is 
interesting to compare the binodal curve in the (r/c? Vp) plane with that of the AOV model 
presented in Ref. ll| (see their Fig. 3). The curve here is significantly more symmetric and 
shifted towards larger values of 77^. 

Finally, we determine the critical point by using the standard cumulant method (we 
essentially follow Ref. 32] )• The order parameter associated with the critical transition can 
be identified as [32 1 

m = 7]^ - {lie) , (14) 

where (rjc) is the average colloid volume fraction. Then, we consider the Binder cumulants 
at coexistence 



M{zp, L) 



U{zp,L) 



(15) 
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FIG. 3: (Color online) Cumulant ratios M and U at coexistence as a function of r/p^*^ for several 
values oi L/Rc- 



The cumulants for different values of L are reported in Fig. [3] as a function of vj'p'^'^ for three 
different values of LjR^. The critical point is given by the intersection point. Both cumulants 
have approximately the same intersection point, which allows us to estimate 



= 129(2), 



/p,crit 



1.321(4). 



At the intersection M and U assume the values 



M„it = 1.226(5), [/ent = 1.57(2), 



(16) 



;i7) 



which are close to those appropriate for the Ising three-dimensional universality class: Merit ~ 
1.239 351 and Ucvn = 1.6036(1) [36^. The small differences we observe are due to field-mixing 



effects 
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381] , which we have not taken into account in the present analysis. At the critical 



point we obtain 



r/c,crit = iVc) = 0.22(1), 



Vp,cnt = (Vp) = 0.93(2). 



Note that these estimates are very different from those of the AOV critical point UJ, ll2| : 



r?e,crit = 0.1340(2), 



Vp,cvit = 0.3562(6). 



(19) 
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^c,crit 
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Ref. 
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0.49 0.21(1) 1.00(5) 
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Ref. 
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Expt. 


Ref. 
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0.11 


0.5 


Expt. 


Ref. 


[42] 


0.92 


0.195 


1.21 


Expt. 


Ref. 
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MC 


Ref. 




0.34 


0.20 


0.29 
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0.19 
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0.18 


0.51 


This work 






0.80 0.22(1) 0.93(2) 



"n 

TABLE II: Experimental (Expt) and Monte Carlo (MC) estimates [4^] of ?7c,crit and r/p^crit- An 



extensive list of results can be found in Sec. 8 of Ref. 



45| 



Both the colloid and the polymer critical volume fractions increase quite significantly. 

The estimate ( fTSi) of the colloid critical volume fraction agrees with all experimental and 
theoretical estimates, see Table Ull for a list of results. On the other hand, the estimate of 



the polymer volume fraction is not consistent with the results of Bolhuis et al. which 
used a much more accurate description of the pair interactions. The discrepancy is probably 
due to our choice of simplified interactions and, in particular, to the fact that polymers and 
colloids can overlap with a relatively small energy penalty. Clearly, quantitative predictions 
require a much more accurate modelling of the interactions among polymers and colloids. 
We also report some experimental results which show, however, large discrepancies. They 



32|: 



are of little use for quantitative comparisons. We finally mention the results of Ref. 
^p'crit — 1-282(2), r^ccrit = 0.150(2), and r^p^crit = 0.328(2). As we already discussed, this 
model is more appropriate for polymers close to the 9 point and, indeed, the estimates of 
the volume fractions at criticality are close to the AOV ones. 

The critical point can also be approximately determined from the behavior of the coexis- 
tence diameter given by 

V/c,diam) '/p,diamj — I 2 ' 2 / ' \^) 

If we consider the intersection of the diameter with a simple interpolation of the coexistence 
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data, we obtain rjc^cnt — 0.21 and ?7p,crit — 0.89. The critical colloid fugacity is consistent 
with estimate ( lT8l) . while the critical polymer fugacity is only slightly underestimated. Size 
corrections, which are not taken into account in this simpler approach, are apparently small 
in the bulk, even at the critical point. As we shall see in the next section, this is no longer 
the case in the presence of the matrix. 



IV. DEMIXING IN THE PRESENCE OF A POROUS MATRIX 

We now study the demixing in the presence of the matrix for a mixture with q = 0.8. We 



take disorder parameters analogous to those used in our previous work [18| : we consider / = 
0.4 and two different sizes of the quenched colloids, Rdis/ Rc = 0.2 and 1.0. Correspondingly, 
Cdis-Rc = 0.070,0.014 in the two cases, respectively. To determine size effects we^ simulate 
systems with L/Rc = 12, 14, 16. We use the same algorithm used in the bulk [18], setting 
A^iter, the number of iterations during which histograms are measured, equal to 40000 A^^- 
Each histogram is averaged over 400 different matrix realizations. Note that system sizes 
are smaller than those used in the AOV case. This is due to the fact that here we are 
dealing with larger polymer volume fractions, hence simulations of the interacting model are 
significantly more time-consuming than AOV simulations of systems of the same size. This 
significantly limits the size of the systems we can simulate. 
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FIG. 5: (Color online) Fluid-fluid binodal curves for / = 0.4, i?dis/^c = 0.2 (left), / = 0.4, 
Rdis/Rc = 1-0 (right). We report the results for L/Rc = 12 and L/Rc = 16 in terms of rjc and r]p. 

A. Demixing curves 



In order to determine the coexistence line Zc = z*{zp) in the {zc, Zp) plane, we analyze the 
colloid and polymer averaged histograms defined in Eq. ( ITOi) using the equal-height method 
(completely equivalent results are obtained using the equal-area method, see Ref. jlSl]): 
coexistence is defined as the value of Zc such that the colloid (or polymer) histogram has two 
peaks of equal height. The histograms at coexistence for R^i^/Rc = 1 are shown in Fig. IHfor 
L/Rc = 14 and several values of r/p'"^. As expected, as r/^'^'^ increases, the two peaks become 
more pronounced, r/cgas (the colloid volume fraction in the colloid-gas phase) decreases, 
while ?7c,iiq (the same quantity in the colloid-liquid phase) increases. As in the bulk case, 
the behavior of the polymer histograms is more peculiar, since the polymer volume fraction 
apparently increases in both phases, except close the the critical point. Again, we expect this 
to be an artifact of the model, due to the simplifications we have introduced. It is interesting 
to stress a second difference with respect to the AOV case. While the AOV histograms are 
strongly asymmetric, with a broad colloid-gas peak and a narrow colloid-liquid peak, in the 
presence of interactions the two peaks are much more symmetric. As a consequence, the 
different methods we used to determine coexistence, the equal-height and the equal-area 
methods (see Ref. ISj for the precise definitions) give fully consistent results, both for the 
colloid and polymer densities at coexistence and for the value z*, for all values of L. 

In Fig. [5] we report our results for the demixing curves in terms of rjc and rip for two 
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FIG. 6: (Color online) Fluid-fluid binodal curves for / = 0.4, Rais/Rc = 0.2 (top) and / = 0.4, 
Rdis/Rc = 1-0 (bottom). We report the results for L/Rc = 12 and L/R^ = 16 in terms of rjcrjp^'^ 
(left) and T]c,r]p (right). 

values of L/Rc. For Rdis/Rc = 0.2 size corrections are small and thus our data allow us 
to determine reliably the infinite- volume coexistence curve. For Rdis/Rc = 1 we observe 
instead some size effects on the colloid-gas branch of the binodal: at a given rjc, the value 
of rip along the binodal increases as L/R^ is increased from 12 to 16. On the scale of the 
figure, the change is small: quantitatively it amounts to an increase of the order of 5-10%. 
In Fig. Owe show the fluid binodals in the reservoir representation, using both r^^''^^ and 77^, 
as we did in the bulk case. For R^is/Rc = 0.2 size corrections appear to be under control, 
except close to the critical point. For R^is/ Rc = 1 corrections are instead larger, especially 
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FIG. 7: (Color online) Estimates of z*R^ as a function of the reservoir polymer volume fraction rf^ 
(left) and of vj^* at coexistence in terms of the reservoir colloid volume fraction rf^. (right). Results 
for L/Rc = 14. 

for the colloid-gas branch, which cannot be reliably determined even far from the critical 
point, where size corrections should be smaller. Clearly, accurate estimates require much 
larger values of L/Rc. 

In Fig. [7] we report z^R^^. at coexistence as a function of rj^ and also rf^* at coexistence in 
terms of the reservoir colloid volume fraction 77^, which represents the volume fraction at the 



given value of Zc-, in the absence of polymers and matrix 46|. Note that, on a logarithmic 
scale, the quantity z^R^^ lies quite precisely on a straight line, indicating that the colloid 
chemical potential at coexistence is well approximated by a linear function in r^^. This 
feature was already observed in the AOV case and seems to be a general feature of this type 
of systems. Fig. [7] differs significantly from that obtained in the AOV case. There, the binodal 
curves showed a significant dependence on /, while here the bulk curve approximately falls 
on top of those corresponding to / = 0.4: they only differ for the position of the critical 
point. This implies that a colloid-liquid bulk phase is almost always in equilibrium with a 
colloid-liquid phase in the matrix, and so does a colloid-gas phase: no capillary condensation 
or evaporation is observed except in a very tiny parameter range, i.e. for those (77^, Vj^) that 
belong to the tiny region between the bulk and matrix binodal curve. 
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B. Critical point 



We wish now to estimate the position of the critical points. This is not an easy task in 
these systems, since the transition belongs to the universality class of the random-field Ising 
model (RFIM) 14j4l7| . Size corrections are large and the efficient cumulant method we used 
for the bulk case does not work. Therefore, we must adopt a different strategy: we follow 
closely Ref. [l5|. First, for each matrix realization we define the averages (Ai'^)vap and (A^^ )iiq 
in the vapor and liquid phases. They are defined as 

(A^')vap = ^ / ^(A^c,int - N,)PiN,)N^, 
-'^vap J 

{N^)u^ = [ e{N, - N,^,,)P{N:)NI (21) 

Here P{N^ is the histogram of N^. at coexistence for the given matrix realization, A^^c^nt 
gives the position of the minimum between the two peaks, 9{x) is Heaviside step function 
[9{x) = 1 for X > and 9{x) = for x < 0], and -ft"vap and K^q are normalization factors. 
Then, for each phase, we define the connected susceptibility jlS] 

X ^ (22) 

where [■] is the average over disorder and V is the volume of the box. Essentially, Xvap and 
Xiiq measure the average squared width of the liquid and vapor peaks for each given sample. 

In Fig. |8] we report the connected susceptibilities for / = 0.4, R^^^/Rc = 0.2,1 and 
L/Rc = 12,14,16. While the results for R^is/Rc = 0.2 are reasonably smooth, the data 
for Rf^is/Rc = 1 are scattered with quite large error bars (we do not even report the data 
at L/Rc = 16, since they would make the figure unreadable). This is probably due to the 
small number of samples we are using. Indeed, in Ref. [i^ it was shown that sample-to- 
sample fluctuations may be quite large and require a particularly large number of different 
matrix realizations to be controlled. In their case, they averaged over 2000-3000 different 
realizations, a number which is significantly larger than ours: we only have 400 different 
samples. In spite of the large errors, the data are, at least qualitatively, in agreement with 
the expected behavior: x(L, rj^) increases with L at fixed rj'p, while, at fixed L, it first increases 
and then decreases as a function of rj^, as also observed in Ref. |15 |. 

General renormalization-group arguments indicate that, close to the critical point, the 



18 



4r 

3.5 



■2.5* 

2- 
1.5- 



1- 



« Rj,«^=0.2L/R=12 
■ Rj, «^=0.2 L« =14 
^ Rj, «^=0.2 L« =16 



i 



I 



1.4 1.45 1.5 1.55 1.6 1.65 



*R,, /R=0.2 L/R=12 

■ dis c c 

# R,. /R=0.2L/R^16 

~ dis c c 



1.4 1.45 1.5 1.55 1.6 1.65 



3.5 



3- 



2.5 



Rj,«^=1.0L«=12 
■ R,,,/R=1.0L/R=14 



1 ' 

i i i i 



9 I I I I I I I I I I I I 

1.5 1.55 1.6 1.65 1.7 1.75 1 



2.5 



1.5 



hi 



i 



R,, /R=1.0L/R=12 

dis c c 

R,. /R=1.0L/R=14 



1.5 1.55 1.6 1.65 1.7 1.75 1. 



FIG. 8: (Color online) Connected susceptibilities x iii the colloid-vapor (left) and in the colloid- 
liquid (right) phases for / = 0.4, Rdis/Rc = 0.2 (top) and Rdis/Rc = 1 (bottom). We report the 
results for L/Rc = 12, 14, 16 as a function of r]p. 



susceptibility satisfies the finite-size scaling Ansatz 



X 



{L/R,)^''F[t{L/R,)^ 



/-I 



(23) 



where 7 and v are critical exponents, F{x) is a scaling function, and t = rj^/rip^^^^^ — 1 
measures the "distance" from the critical point. This scaling Ansatz is, strictly speaking, 
valid in the limit L — )■ 00, t — )■ at fixed t^L/RcY^'^ and neglects (analytic and nonanalytic) 
scaling corrections. 



As discussed at length in Refs. 



19| . the critical transition belongs to the RFIM 
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V 



Newman & Barkema [47] 1.02(6) 0.15(7) 
Hartmann & Young [48] 1.32(7) 0.50(3) 
Middleton k Fisher [4^ 1.37(9) 0.51(3) 
Vink et al. [15] 1.1 0.13 

Fernandez et al. [50] 0.90(15) 0.531(40) 



TABLE III: Estimates of the critical exponents u and rj for the RFIM universahty class. The 
exponent 7 is related to v and r] hy ^ /v = 2 — r]. 

universality class. Some numerical estimates of the RFIM critical exponents are reported 
in Table IIIII It is evident that there is no general consensus on the estimates and thus we 
tried all different possibilities. Let us first consider the data for R^i^/Rc = 0.2, which are 
the most precise. We determine the critical polymer fugacity Zp,crit, or equivalently ?7pcrit5 
by requiring the estimates of xL~^^'^ to collapse onto a single curve as a function of tL^^'^, 
fixing the exponents to the RFIM values. If we use the numerical values of the exponents 



employed in Ref. 15[, we obtain a poor collapse of the data for any choice of the critical 



polymer fugacity. Apparently, the best collapse is obtained by taking the latest numerical 



estimates of Ref. [50|. The rescaled data are reported in Fig. IH] and allow us to estimate 
the scaling function F{x) appearing in the scaling Ansatz (123|) . This function is universal 
apart from two rescalings: if Fi{x) and F2{x) are determined for two different transitions 
belonging to the same universality class, then one should have Fi{x) = (6x), for suitable, 
nonuniversal constants a and b. We can thus compare the curves appearing in Fig. [9] with 



the analogous ones reported in Ref. 15[. Shapes are quite similar, although in our case the 



peak of the scaling functions apparently occurs for x ~ 0, while in Ref. [15| the maximum 
occurs for x well below zero. The origin of this difference is unclear: it might be due to 
the different choices for the RFIM critical exponents or to the neglected scaling corrections. 
The results for the critical parameters are reported in Table IIVI It is interesting to compare 
these results with those which would be obtained by a more naive analysis of the diameters. 
If we compute the intersection of the diameter line with an interpolation of the binodal data 
for L/Rc = 14 we would obtain T^p^crit = 

1.02, 

Vc,cvit — 0.20. The critical colloid volume 
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FIG. 9: (Color online) Rescaled connected susceptibilities (L/Rc) ^^^x ^ ^ function oft {L/R, 



where t = Vp/^p,crit " 1' fo^' 
Rdis/Rc = 1 (bottom). We use 



12, 14, 16. Data for / 
50t] ^/zy = 1.5, ly = 0.9. 



0.4, Rdis/Rc = 0.2 (top) and 



fraction agrees within errors with that reported in Table llVt while rjp^crit is slightly (6%) 
underestimated . 

It is not easy to extend the analysis to the case Rdis/Rc = 1, given the very large fluctu- 
ations in the data. A reasonable collapse, see Fig. [9l is obtained for the parameter choices 
given in Table IIVI However, errors cannot be reliably estimated. Note that size corrections 
are quite significant for this size of the quenched colloids. For instance, half of the points 
for which we observe bimodal histograms belong to the homogeneous phase in the infinite- 
volume limit (compare the estimate of 1]^ ^^.^^ with the data reported in the right-bottom panel 
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,crit ^c,crit ^p,crit ^c,crit 



bulk 1.321(4) 0.575(1) 0.93(2) 0.22(1) 

0.4 0.2 1.49(6) 0.602(4) 1.08(2) 0.21(1) 
1.0 1.67 0.62 1.17 0.22 



TABLE IV: Estimates of the critical point position. 
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FIG. 10: (Color online) Binodal curves for the interacting-polymer model (solid symbols) and for 
the AOV model (lines). Empty symbols give the critical points. 

of Fig. ED. 

V. CONCLUSIONS 

In this paper we consider a simplified CG model for polymer-colloid mixtures. Polymers 
are represented by point particles interacting by means of pair potentials. Models of this 
type have been shown to be quantitatively predictive as long as the polymer density is below 
the overlap density, i.e. for rjp < 1. Hence, they can be used in the study of mixtures in the 
colloid regime {q ^ 1), in which fluid-fluid phase separation occurs for rjp < 1. We study 
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polymer-colloid segregation in the bulk and in the presence of a porous matrix for / = 0.4 
and for two different values of Rdis/Rc, 0.2 and 1. In Fig. [10] we summarize our results for 
the coexistence curves and compare them with the AOV ones. The effect of the matrix is 
similar in the AOV and in the present interacting model. Disorder moves — but the effect is 
not large — the coexistence curve towards larger polymer densities. In the colloid-gas phase 
(for rjc = 0.1, say), the observed change is quantitatively similar in the two models. The 
volume fraction rj* at coexistence increases by 10-20% both in the interacting and in the 
AOV case as / increases from / = (bulk) to / = 0.4. The behavior in the colloid-liquid 
phase is instead quantitatively slightly different. In the AOV case, the /-dependence of the 
coexistence curve is very small (all curves coincide within errors), while in the interacting 
case the behavior is analogous to that observed in the colloid-gas phase: i]* at coexistence 
increases by 10-20%. It is important to stress that this different behavior may be an artifact 
of the model, since, as we have already observed several times, binodals cannot be trusted 
for ?7c > 0.35. 

While binodals do not change significantly with disorder in both cases, the critical point 
position changes significantly. Moreover, the behavior is quite different in the AOV and in 
the interacting-polymer model. In the AOV case the critical colloid volume fraction increases 
significantly as / changes from zero (bulk case) to 0.4 and 0.7, while r^p^crit is approxim ately 
constant or slightly decreases: for / = 0, r^ccrit = 0.1340(2) and rjp^cnt = 0.3562(6) ll|, Il2[ |. 
while for / = 0.4 r^ccrit ~ 0.19 and ?7p,crit varies between 0.27 and 0.34 depending on R^is/Rc 
[isl . In the interacting case instead, the colloid recent is approximately independent of /, 
while rjp^crit changes — it increases — with both / and R^is/Rc- Note that a similar behavior is 
observed for the q dependence (at least, in the colloid regime) of ?7c,crit and rjp cnt in the bulk 
^. In the AOV is essentially independent of g, while //ccrit varies significantly 

with q. In the interacting-polymer case the opposite occurs: r^p^crit varies with q, while r^ccrit 
stays approximately constant. 

In the AOV case we found lisj] that colloids could undergo capillary condensation: for 
some values of the parameters, a colloid-gas phase in the bulk was in chemical equilibrium 
with a colloid-liquid phase in the matrix. In the present model no such phenomenon occurs 
unless one carefully tunes the system parameters. 

The results we have presented rely on simulations of a very simplified model: first, we use a 
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CG model of the polymers in which each polymer is replaced by a monoatomic molecule; sec- 
ond, we use simplified potentials which allow us to reproduce correctly the thermodynamics 
in the low-density limit, but not the intermolecular distribution functions. An improvement 
of the model is clearly necessary if one wishes to obtain quantitatively accurate predictions. 
A multiblob model |25| in which polymers are represented as polyatomic molecules is neces- 
sary if one wishes to consider larger polymer-to-colloid ratios or to obtain accurate results for 
small values of Rdis/Rg- In the latter case, one would expect that an accurate modelling re- 
quires CG multiblob models in which the blob size is smaller than -Rdis- Studies of the phase 
behavior of mixtures in random porous matrices with these improved models are not feasible 
with present computer resources. However, a less CPU-demanding task would be the study 
of the phase behavior in the presence of regular arrays of obstacles (the quenched colloids 
could belong to the sites of a regular lattice). Although less interesting from an experimental 
point of view (silica gels have a random distribution of pores) these type of systems could 
be more carefully studied, without relying on many different approximations. They can thus 
provide a theoretical laboratory, where one can understand the role of quenched obstacles 
on the phase behavior of these systems and thus develop theories that can then be extended 
to the more difficult random case. Work in this direction is in progress. 
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